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Abstract. We discuss tiie many-body piiysics of an ensemble of Rydberg dressed 
atoms with van der Waals dipole-dipole interactions in a one-dimensional lattice. 
Using a strong coupling expansion and numerical density-matrix renormalisation group 
simulations, we calculate the many-body phase diagram. A devil's staircase structure 
emerges with Mott-insulating phases at any rational filling fraction. Closed analytic 
expressions are given for the phase boundaries in second order of the tunnelling 
amplitude and shown to agree very well with the numerical results. The transition 
point where the incompressible phases melt due to the kinetic energy term depends 
strongly on the denominator of the filling fraction and varies over many orders of 
magnitude between diff"erent phases. 
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1. Introduction 

Recently many-body systems with non-local, power-law interactions gained considerable 
interest as the non-local coupling can give rise to quantum phases that do not exist for 
point-like interactions [Ij. For example in one spatial dimension repulsive dipole-dipole 
or van der Waals interactions of atoms can lead to a variety of crystalline ground state 
phases with less than unity filling in the presence of a commensurable periodic lattice 
[21 El HI [3| . For very strong lattice confinements the filling fraction forms a so-called 
complete devils staircase as function of the chemical potential /x, i.e. every rational 
filling between and 1 is stable for a finite interval of /x and these intervals form a 
dense staircase [^. Power-law interactions arise e.g. for dipolar atoms or molecules 
[1] . A very interesting alternative approach are Rydberg gases for which there has been 
considerable experimental progress in the recent years [3 131 ^ [IDl CH [El [E] . While most 
previous experiments implemented schemes where Rydberg excitations were created by 
continuous near-resonant laser driving, alternative proposals have been put forward to 
use far-detuned light fields to admix a small component of a Rydberg state to the ground 
state of atoms [21 [HI [15]. The potential advantage of this "Rydberg-dressing" is that the 
interaction strength can be tuned to a certain extend and is reduced to an energy scale 
where the center-of-mass motion of atoms can become relevant. Furthermore continuous 
driving by a near-resonant laser does in general not conserve the number of Rydberg 
excitations, while in the case of Rydberg-dressing the number of atoms are to good 
approximation a conserved quantity. While the non-local repulsive interaction favors 
crystalline structures with non-unitary filling in lattice gases, hopping between lattice 
sites induced by tunnelling of the atoms can lead to a melting of the crystals [El [E] • 

In the present paper we analyse the melting process induced by hopping of 
Rydberg-dressed atoms. In particular we derive the phase boundaries of stable 
crystalline structures as a function of the hopping amplitude J both analytically and 
with exact numerical simulations. To this end we employ a second order strong 
coupling expansion on the one hand and numerical density-matrix renormalisation group 
(DMRG) simulations, adapted to long-range interactions, on the other. Both show 
excellent agreement up to second order in J . A corresponding strong-coupling analysis 
of the phase diagram has been discussed before for the case of dipolar interactions 
in [6], however no explicit expressions were given. We here derive analytic expressions 
which in the case of fast decaying potentials, such as the van der Waals potential, 
can in very good approximation be written in a compact closed form. Furthermore 
we provide for the first time a comparison to numerical data from DMRG simulations 
adapted to long-range interactions. 

As a starting point we consider the generalised Bose-Hubbard model in one 
dimension 

^ . ' ^ — . ' 

^0 
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For convenience we have set = a = 1, where a is the lattice constant. hi^h\ are the 
annihilation and creation operators at site and J is the hopping amplitude between 
neighbouring lattice sites describing the tunnelling of atoms between adjacent lattice 
sites [18j. Most parts of our discussion are valid for any convex interaction potential, i.e 

V{r + l) + V{r-l) >2V{r). (2) 

However to be specific we discuss in the following power-law potentials of the form 

with /3 G N\{1} and (7/3 > being the interaction coefficients. /5 = 6 corresponds to 
a van der Waals type interaction which results from virtual (i.e. off-resonant) long- 
wavelength electromagnetic transitions in the manifold of Rydberg states. If there is an 
accidental or engineered resonance, a so called Forster resonance, also the case /5 = 3 
can be of relevance. The model interaction potential ^ diverges for r ^ and thus 
two particles cannot sit on the same lattice site. To very good approximation the same 
is true if one considers the exact interaction potential between Rydberg atoms which 
deviates from the above power law for small distances (see e.g. [191 [20]). f^ct for the 
Rydberg-dressing scheme one finds an effective two-body interaction potential between 
atoms at positions and [21 [Ml [15] 

'^^'-''^= \r.-r,\^ + Rr ^^.^ r = h-r.|. (4) 

Here Cq is the Rydberg interaction coefficient and Rc = {Cq/{2A)Y^^ describes a 
characteristic length scale of the interaction potential, where A 1^2 1 is an effective 
detuning of the off-resonant laser driving, ry ^ f]^/A^ <C 1 accounts for the small 
admixture of the Rydberg state. Considering e.g. ^''Rb with a Rydberg state of principle 
quantum number around n = 60 and detunings of a few tens of MHz one finds a value 
of Rc on the order of one to a few fim. In a recent experiment atoms were loaded into 
optical lattices and excited to Rydberg states with n = 55 — 80 [2JJ. In this experiment 
the lattice spacing was tunable between 0.5//m and 13/xm. Thus the cut-off length Rc 
may become smaller than the average distance between particles given by the lattice 
constant devided by the average occupation number per lattice site. In this case its 
effect can be disregarded and this is the parameter regime we are considering here. For 
that reason 6^, are assumed to describe hard-core-bosons with commutation relations 

Kb]] =[kh] =[^U]]=0 tj^j, (5) 
{bM={blbl} = 0, (6) 

{h:bi} = ^- (7) 

We note that although we discuss here the particular case of the integer value /3 = 6, 
the qualitative structure of the phase diagram holds true for any positive value of /3. 
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Figure 1. Ground states for filling fractions q = \^q= \ and g = | at J = 0. 



2. Ground state for J = — classical limit 

a.marueg@web.de Let us first consider the ground state of Hamiltonian ([T]) in the limit 
of vanishing hopping, i.e. J ^ 0. In this limit the Hamiltonian is diagonal in the 
number basis, where each lattice site contains exactly zero or exactly one particle, 
which corresponds to a classical situation. Consequently the energy can be minimised 
just by finding the configuration of classical particles that gives the smallest interaction 
energy. As the interaction potential is convex (see equation Q), the minimum energy 
for commensurable filling fractions q = ^ < where m, n G N are relatively prime, is 
attained by a regular pattern with unit cells of size n [22] . For a given chemical potential 
the ground state is such a phase with rational q. For these values of q the corresponding 
phase is incompressible, i.e., particle as well as hole excitations require finite energy, i.e. 
/x+ > where 



= ±{E^{q)-E{q)). (8) 

p=q±0 

Here E{q) is the ground state energy for filling fraction q. E'^{q) is the corresponding 
energy where one particle has been added {E^[q)) or respectively removed {E~{q)). p 
is the average number of particles per lattice site. Because of particle- hole symmetry, 

fiiq) = n-li{l-q), (9) 

where Q is the energy per particle in the completely filled lattice is the Riemann zeta 
function ({/3) in the case of ([s])), it is sufficient to consider only filling fractions ? < |. 

The phase-diagram in zeroth order of the hopping J has been discussed some time 
ago by Bak and Bruinsma [23] . Here, we will brieffy review their results and we will 
make use of their notation: denotes the position of the i-th particle. Xf = X^ — X^ 
describes the distance between the p-th and the i-th particle where the particles are 
numbered from the left to the right. In the ground state of Hubbard's solution [22] 
shows that all possible distances are given by 

Xf = rp or Xf = rp + 1, (10) 

where < | < + 1. If | G N then Xf = = | and all particles are equivalent. 
Examples of ground states are shown in figure [TJ Due to the convexity of the interaction. 
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Figure 2. Defects for q = = | and ^ = |. The boxes mark the extension of the 
defect. Outside of the boxes the configuration of particles is exactly the same as in the 
commensurable case, see figure [l] 



all nearest-neighbour separations are maximally close to the average separation, i.e. they 
are either [^J or [^] . 

The chemical potential for J = is given by 

^i^2\q) = ± {{q^\Ho\q^) - {q\Ho\q)) . (11) 

\q) denotes the ground state of Hq for filling fractions g, whilst denotes the 

corresponding ground states where one particle has been added respectively 

removed \q~)- Although both expectation values in (11) are infinite in the 



thermodynamic limit, their difference is finite, and can be calculated by summing the 
change in interaction energy particle by particle. The same is true for the first and second 
order corrections in the hopping amplitude J discussed below. For filling fractions q = ^ 
there will be n defects for because it is energetically favourable to break up an 

extra particle (or a hole) into n defects each with fractional charge Some of these 
defects are displayed in figure [2] In the thermodynamical limit they will be separated 
by arbitrarily large distances so that we can assume without approximation that they 
do not interact with each other at all. Note that is not uniquely defined, but the 
position of defects is arbitrary, as long as they are well separated. For filling fractions 
of the form q = - the chemical potential is given by [23] 



(12) 



We can go a step further and evaluate (12) for the power-law potential ^ analytically. 



n/^(/5 - 1)! 



n(/3-i)(/3-2)! 



n 



(13) 



^^-^ and its l-th derivative "i/^^^z). To describe more 



using the digamma function ^^(2:) — ^^^^ 
general filling fractions q = ^ with m 7^ 1 as well, we only have to replace ^ + 1 
for /x^^ in (12) if N. Plotting the filling fraction as a function of the chemical 
potential gives rise to a complete devil's staircase (see figure [3]). Every rational filling 
fraction between and 1 has a finite stable range with respect to the chemical potential 
and the functional dependence consists of a dense set of steps between these stable 
plateaus, a so-called devils staircase. The union of these intervals covers the full range 
of /X. It is clearly seen that filling fractions with small denominator n are the most stable 
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a) /''°V\ b) ^''"VAs 

Figure 3. Devil's staircase: The filling fraction q in the ground state J = is plotted 
versus the chemical potential ji^^^ for a van der Waals interaction potential /3 = 6 (a) 
and a Forster resonance /3 = 3 (b) in units of the level shift = C^/a^ . The insets 
to the upper left emphasise the repeating features of the devil's staircase. The insets 
at the right bottom show this equation of state in a double logarithmic plot. The red 
line in (a) and the green line in (b) correspond to q = {ji^^^ / C ^ . 



as can be seen directly from equation (13). These intervals of stability depend crucially 
on the exponent (5: for large /5, the interval of the stable phase q = \ overgrows by far 
the size of all the other phases. Note that for n big enough (n=F l)/n is almost constant 
and hence the mean chemical potential ji^^^ = (/x^^ — /x^^)/2 scales as . Thus below 
half filling the fraction q scales as {jji^^yC{3)~^ as indicated in the inset of figure [sj 

As mentioned in the introduction the true interaction potential between Rydberg- 
dressed atoms becomes fiat below a certain distance. As a consequence the atoms can 
be treated as hard-core bosons only for sufiiciently small energies, i.e. sufiiciently small 
chemical potentials. The finite cut-off Rc will lead to a qualitative change of the phase 
diagram in parameter regions where the separation between atoms becomes smaller than 
Rc. This is however only the case if the filling q becomes large and in particular exceeds 
the ratio Rc/a. 



3. Perturbation theory in J 

We are now interested in the melting of a crystal phase with increasing hopping rate J . 
Therefore we perform perturbation theory up to second order in J . Similar calculations 
have been done by Burnell et al. for dipolar interactions ^ but no compact 

analytic expression was given. 

Let us start with first order processes in J . Then, the chemical potential reads 
IJi±{q) = ± (^{q^\Hj\q^) — {q\Hj\q)^. Hopping of any single particle in state \q) will 
not contribute to any energy correction, because the resulting state has no overlap with 
the J = ground state. The same is true for hopping of any but the 2n particles of the 
state \q^) which sit at the left and right borders of any of the n defects. Hopping of 
these particles will lead to a hopping of the respective defect by n-sites, see figure [4j As 
the states with localised defects are degenerated with respect to /fo, we have to apply 
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Figure 4. Hopping of defects for g = ^^q = | and ^ = §• The black arrow indicates 
the hopping of one particle and the red one corresponds to the resulting hopping of 
the defect. The blue ones are the alternative hopping-possibility, which would result 
in a hopping of the defect in the opposite direction. 
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Figure 5. All possible second order processes for 1^^+) (a) and \q~) (b) at filling 
q = here ^ = | • (i) generates an effective chemical potential for single particles that 
depends on the distance to the defect, i.e., the defect polarises the background, (ii) 
refers to a virtual deformation of the defect, (iii) corresponds to an effective hopping 
of the defects over two unit cells. 



degenerate perturbation theory. Therefore, we look for a basis of these states in which 
Hj is diagonal. This is given by states where all the defects are delocalised over many 
lattice sites (yet well separated from each other) with some quasi-momentum k. The 
energy is minimised for a state with quasi-momentum zero. The resulting first order 
correction to the phase boundary thus takes the simple form 

/i2^(g) = T2nJ. (14) 

For small but nonzero J gaps open between the incompressible phases of rational 
filling, because the kinetic energy gained by delocalisation of defects favours a finite 
density of defects. At the tip of the phase J becomes large enough such that the 
creation of particle and hole defects becomes favourable even at commensurate filling 
fractions and the crystalline phase melts. This is completely analogous to the ordinary 
Hubbard model. However, for any arbitrarily small but finite value of J almost all 
phases with rational filling become unstable, leaving only a finite number (those with 
the smallest denominator n) and destroying the devils staircase.a.marueg@web.de 

In second order perturbation the energy corrections are given by = 
(^l-^^I^X^I^^I^) ^ where |^^) denotes the ground state of and {|$)} is a complete 

set of orthonormal states. In the following we consider only filling fractions of the form 
g = ^. In principle it is posa.marueg@web.desible to extend the calculations to other 
fractions. However, finding a general formula which describes all energy differences in 
the denominator for all hopping processes is a far more tricky task. In figure [5] all 
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possible processes for = are shown. As has been seen in first order, the states 
with localised defects are degenerate with respect to ^o- In degenerate second order 
perturbation theory the energy correction is given by {q^\HjPHj\q^) /{E^^J — E^'')^ 
where P projects out the subspace of ground states of ^o- The state \q^) has to be an 
eigenstate of Hj restricted to the degenerate subspace. Therefore, intermediate states 
1$) with E^"^ = do not contribute to the energy corrections. For states \q) only 
process (i) in figure [i] contribute. For the commensurate state \q) one finds for the 
second order correction to the energy 

where is the number of particles in the lattice. AE^^\q) = E^^\q) — E^\q) is the 
energy difference between the ground state \q) and the same state but with one particle 
hopped by one site. With the notation introduced above it can be evaluated to 

CO 

AE^°\q) = J2 {^Virp) - V{r, - 1) - V{r, + 1)) (16) 



2C{/3) {-ir 



n J \ n 



n^{/3-l)\ 

In order to obtain the (finite) chemical potential we have to calculate also the energy 
corrections E^'^\q^) for the states with one extra particle or hole \q^) in second order 
perturbation theory. To this end let us have a closer look at all possible hopping 
processes displayed in figure [5} As has been done in first order, we will concentrate 
on a single defect multiplying the resulting contribution with its number n. Process 
(zz) describes a virtual deformation of the defect while process {in) corresponds to an 
effective hopping of the defects over two unit cells. Both processes have only to be 
taken into account for one particle of the defect and one particle on its right side and 
therefore contribute only to a finite energy value. In the thermodynamical limit, process 
(z) takes place for an infinite number of particles. For this reason, its contributions to 
the energy corrections will diverge. However, for a decaying potential the contribution 
to the energy correction of a particle far away from the defect will be the same as the 
contribution of a particle of the commensurate state \q). Subtracting now the energy 
corrections for states \q) and these contributions of process (z) will mostly cancel. 
We find 

where AE^^j^ denotes the energy difference between the ground state \q^) and the same 
state, where the j-th particle right of the defect moves right (second subscript +) or 
left (second subscript — ). For the meaning of position see figure [5j We have 

assumed here that the defect is in the centre of the system and thus the summation 
limit is = l/2[(A±l)g — 2]. As \q^) has to be an eigenstate of Hj in the degenerate 
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subspace the defects are completely delocalised. As in first order it turns out that the 
ground state has quasi-momentum k = 0. Then, the second order correction to the 
chemical potential can be written as 

1 1 2 
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A<5 
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1 - liSn,2 ± Sn,2) 
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(18) 

(19) 
(20) 



Part (18) takes into account all processes of type (i). Because of symmetry it is only 
necessary to sum over particles on the right side of the defect. The last term in this 
line is due to the fact that there are less particles in Iq"^) taking part in process (i) then 
there are particles in state \q). Part (19) describes process (ii) and part (20) process 
(in). Calculating all the energy differences analogous to the case of AE^^\ the chemical 
potential for filling fractions q = - can finally be written in second order in J as 
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(21) 



where the terms Sf{n) are given by 
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(22) 



If the power (5 of the interaction potential is sufficiently large and if the effect of the cut- 
off length Rc can be disregarded, i.e. if Rc < a/q^ the entire sum S^{n) contributes 
only very little. For the case of a van der Waals potential, i.e /5 = 6, we have numerically 
verified that the sum contributes less than a few per cent to the final result. Thus we 
can to a good approximation ignore the sum, which gives for the chemical potential of 
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the q = 1/n phases in a system with van der Waals interactions 

(2) f 7560n6(2T -) 

^± Cp\ l67r3-63[*(^)(^) + *(5)(^)] 

^ 7560n6(l-i(5„,2±M) 



87r3 + 63[*(5)(^^) - *(^)(^) - 
+ 2 / — V (23) 

n6 (n+l)6 J 



4. Exact numerical calculation 



Complementary to the perturbation analysis we perform numerical calculations for the 
case of a van der Waals potential. For this we employ the density-matrix renormalisation 
group algorithm [24] . It is a variational technique that minimises the energy of the full 
Hamiltonian Hq + Hj using a matrix product state (MPS) ansatz [25j. Although limited 
in the amount of entanglement along the lattice that they can capture, MPS are known 
[26] to approach the true ground state of one-dimensional systems quickly with growing 
matrix dimension, if interactions are of finite range. Although the interactions decay 
only polynomially in our model, we find that we can safely cut-off the interaction at 
a finite distance of r lattice sites, and the ground state energy will be for all practical 
means independent of r as long as we restrict the filling fractions to denominators n that 
are small compared to r. In order to make our DMRG implementation more efficient 
we group the interaction terms as originally introduced for computations in momentum 
space [27j. We also take advantage of the particle number conservation, which allows us 
to use MPS with the proper symmetry and fix the total particle number a priory, which 
comes in handy to calculate the phase diagram. 

To calculate the phase boundaries for a given filling fraction we make direct use of 
([s]) by computing E{q) and E^{q) for a finite system[f| The minimum system size L to 
capture the physics in the thermodynamic limit for a given filling fraction q = ^ can be 
estimated to be 2n^, because it must fit n defects of size about n separated from each 
other and the boundaries by at least one unit cell of size n. For not too large n we can 
perform a infinite size extrapolation by employing different system sizes and assuming a 
1/L scaling of the finite size error. In figure [6] we plot the resulting phase diagram for a 
van der Waals potential. The dashed lines display results from perturbation theory 
including up to the second order and the continuous lines are DMRG results. As 
expected the agreement is very good for small J/Aq. Perturbation theory however 
fails close to the critical points corresponding to the tips of the lobes of incompressible 
phases. Also using DMRG the exact position of the critical point is hard to compute due 
to strong finite size effects. Accordingly we have not attempted to extract these values. 

I We present results for open boundary conditions, i.e. no tunnelling or interaction between the left 
and the right end of the system, here. Periodic boundary conditions are also possible, but the lack of 
boundary effects does not make up for the higher numerical cost [28 which limited our computations 
to smaller system sizes. 
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Figure 6. Phase diagram for a van der Waals interaction potential. The dashed hnes 
correspond to perturbation theory up to second order in J/Ae, with Ae = Ce/a^, 
continuous hnes are calculated with the DMRG method using MPSs of dimension 32. 
Coloured lines are infinite size extrapolations for L > 24 (n = 2, 3 only), L > 60, and 
L > 120 lattice sites. No infinite size extrapolation is applied for n = 6 because we 
decided not to make the effort to calculate large enough system sizes. For all cases 
interactions over distances larger than r = 7 lattice sites, a) Double logarithmic plot. 
The finite size results are shown in grey for (7 = | and (7 = | for illustration, b) Linear 



plot. Note the vast difference in size for different n (13), which is much more drastic 
than in the /3 = 3 case [6 . 



The continuous lines in the figure end at arbitrary values, while the phase boundary 
ends where these lines make contact. 

As can be seen from fig jH] the melting points of the different incompressible phases 
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differ by many orders of magnitude in the normalized detuning J/Aq. This raises the 
question if the corresponding timescales are compatible with the finite lifetime of the 
dressed Rydberg atoms. It should be noted however that the relevant time scale is 
determined by Aq = C^/a^ and thus depends on the lattice constant a. As mentioned 
in the introduction the interaction potential is only valid beyond a cut-off length 
which sets a constraint on a. On the other hand the properties of the incompressible 
phases with filling q are only determined by the tails of the interaction potential at 
distances a/q and thus for our model to hold it is sufficient that a > Rcq. (For the 
same reason the smearing out of the interaction potential due to the convolution with 
the Wannier functions has little effect.) This means Aq and thus the relevant time scale 
can be modified by adjusting the lattice spacing according to a ^ Rcq resulting in a 
scaling Aq ^ q~^. For example taking a = Rcq one finds that the melting points of 
the q = 1/2 and q = 1/6 phases are J1/2 ^ 6 x Cq/R^ and Ji/q ^ 0.5 x Cq/R^. Since 
Cc/Rc ^ f^^/A, where Q and A are the effective Rabi- frequencies and detunings of the 
Rydberg dressing, the common energy scale of the meltig points can be tuned and can 
be in the kHz to MHz range. Thus melting can be observed also for small fillings well 
within the lifetime of Rydberg dressed states. 

5. Summary 

We have calculated the ji — J phase diagram for Rydberg dressed atoms in a deep 
one-dimensional lattice potential. For vanishing hopping strength, J = 0, the 
chemical potential ji forms a complete devils staircase as function of the filling fraction 
corresponding to stable crystalline phases at any rational filling. Following a similar 
analysis of reference [23j we derived a closed analytic expression for any power-law 
interaction potential. We then analysed the melting of the Rydberg crystals with 
increasing hopping. To this end we performed a second order strong-coupling expansion 
and found excellent agreement to exact numerical simulations based on DMRG adapted 
to long-range interactions. For the case of stronger localised interactions, such as a van 
der Waals coupling, a compact analytic expression for the phase boundaries of phases 
with filling q = 1/n was found. 
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